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ABSTRACT 


An iterative procedure of Wiener filtering for tfca 
restoration of images proposed in (2/7) is studied in 
detail for the restoration of images degraded by imaging 
systems and further corrupted by additive and multipli- 
cative noise. It is further extended to include genera- 
lised Wiener filtering and suboptima 1 filtering in other 
than frequency domain in order to reduce the number of 
conputations. Simulation results support the effective- 
ness of the iteration strategy and shows that the mean 
square error converges in few iterations (typically three 
to four) and good restorations are obtained. The iterative 
procedure overcomes some of the shortcomings associated 
with Wiener filtering. 
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CHAPTER 1 


IMAGE RESTORATION 


Digital restoration o£ pictures deals with 
images which &xv degraded by imaging systems and further 
contaminated by extraneous noise. In some cases, the de- 
gradation i# a ‘point degradation* process where only the 
gray levels of the individual picture points are affected 
without introducing blur. Other types which involve blur 
are known as ‘spatial degradations 1 . Still further degra- 
dation arises from tenporal or chromatic effect. Here we 
will consider only the first two processes which corrupt 
a picture. These occur in various applications. For 
exanple, in astronomy and remote sensing, the pictures 
obtained are degraded by atmospheric turbulence, ' 
aberrations of the optical system and relative motion 
between the object and the camera. For medical radiographic 
images, both the resolution and contrast are poor due to 
the nature of the X-ray imaging systems. Electron micrographs 
are often degraded by the spherical aberration of the electron, 
lens. One also encounters non-linear imaging systems. However 
restoration in non-linear models requires enormous corrputa— 
tional effort and has not achieved any significant progress. 

An a result the assumption of linearity is ..usually made, 
linear imaging systems may be shift variant/invariant. Further 
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the image can be contaminated by both additive and multi- 
plica tive, signal dependent and independent noise during 
detection and recording process. 


In the following, we will consider the presence 
of both additive and multiplicative noise. 

IMAGE FORMATION AND RECORDING, MODEL : 


n 1 (x,y> 


f (x,y) 



Imaging system. 


g(x,y) 


|n 2 (x,y) 

Detection and recording system. 


Fig. 1.1 Continuous image formation and recording model 

Fig. 1.1 shows a comp rehen sive model ( 2 ) of the 
imaging and recording system for continuous signals (The 
discrete version of this model is quite straightforward and 
will be taken up in Chapter 31. The following notations are 
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used in this discussion. 

£(x,y) = Gray level function of the original, otoj ect. 

* 

k(x,y,x*,y*) « Point spread function/lnpxilse response of 

the imaging system i.e. the output of the imaging system at . 
(x,y) corresponding to an inpulse input described by f(x#y) = 
6(x-tt, y-0) will be given by h(x,y/**fJ)# The PSP (Point- 
spread function) is, in general, dep endent on the position. 
(ct,,0) of the inpulise in the original picture. 

SCO * Detector response function. In the following 
discussion, we have assumed it to be identity transformation!. 

n^ (x,y) =* Multiplicative noise-component. 

n 2 (x,y) « Additive noise -conponent. 

The final output g(x,y) is related to the input 

object by 

g(x,y) = //h(x,y,x , ,y , )f(x*>y f ) dx* dy' 

+ n^Cx/y) // h(x,y,x , ,Y‘)f (X* ,y* )dX , dy , -Ha 2 (x,y) . 

.* a.D 

If> however, the degrading system is space-invariant, then the 
PSP. takes the form h(x-x*, y-y 1 ) and accordingly. r. 
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gC*,y) *J‘/h(x-n , ,Y-y«)f(x* # yMdm*dy*+n 1 (x # y)// h(*-x*,y-y‘) 
f(X* ,y’ )dJ*'dy* + n 2 (x # y) •• (1*2) 

The objective of the estimation procedure is to 
make as good an estimate of the original picture as possible 
given the model (Pig. 1.1) and the degraded picture. Any 
such estimation procedure requires some knowledge concerning 
both the degradation function and the noise. In some cases, 
the physical phenomenon governing the degradation can be used 
to determine the PSP. In other situations, the PSP may be 
found out from the degraded output itself with a test-picture 
at the input ( 4 ) . Apart from the PSF, also required are 
the statistical properties of the noise and the knowledge 
about how it is correlated with the picture. The most common 
assumption about noise is that it is white i.e. its spectral 
density is constant, over all frequencies and it is uncorrelated 
with the picture. The concept of white noise is a mathematical 
abstraction and though not always admissible, it is a reasonable 
assumption provided the noise bandwidth is much larger than the. 
picture bandwidth. 

In what follows, we will start with the restoration 
of one-dimensional signals. Modifications required for two- 
dimensional applications are straightforward and will be explained 
later on. 



CHAPTER 2 


CONTINUOUS ITERATIVE WIENER FILTER 

The restoration of images in a continuous model 
will be considered in this chapter. We will assume the 
degrading system to be deterministic and space invariant* 
Also we will restrict ourselves to the restoration of one- 



Fig. 2,1 Imaging and recording model for one-dimensional 
signals with space-invariant degrading system* ' 

Fig, 2,1 shows the corresponding imaging model with h(x) 
denoting the impulse response of the space— invariant degrading 
system. The output g(x) is given by 

gCx) * / h(x-x* )f (x')dx* + n x (x)<f h(x-x*) f (x* )dx* 

+ n 2 (x) ,, (2.1) 

Fourier-trans forming both sides, 

<3 ( j w) * H ( j.w) . F ( j,w> + N x ( j w) * (H ( j w) ,F (j ; w) ) 
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where G(jw)> H(jw) # F(jw)> % (jjw) and N 2 (jw) are the Fourier- 
trans forms of g(x), h(x) t fix), n-^(x) and (x) respectively 
and ' * * denotes convolution. 

In the absence of noise/ we have 

G(jw) «* H(jw) .P(jjw) 
or equivalently, 

FCjw) » G(j,w)/H(jw) .. (2.3) 

which implies that f(x) can be restored by multiplying the 
Fourier transform of the degraded output by the function 
H ( jw) and then inverse Fourier transforming i.e. the filter 
function is hnown as ‘Inverse filter*. 

However, a number of problems arise when one atterrpts to make 
practical use of it ( 4 ) . There may be points on the w— axis 
where H(jw) is zero and consequently G(jw) is also zero # 
resulting in indeterminate ratios. In the presence of noise, 
the zeros of G(jw) do not coincide with those of H(jw) . Aa 
a result/ the division by H(jw) would lead to very large 
values in the vicinity of the zeros of H(jw). In fact, when 
noise is present, we have 

Gfjw) *H(jw),P(jw) + % (jw) * <H(jw) .F( j ; w) ) + N 2 Uiw) 



or. 


£P4 

H(jw) 


(Nj^Cjw)* (H(jw) .P(jw) ) ) 

+ irrj95 


N 2 (jw) 
+ H(jw) 


.. (2.4) 


When H(jw) is very small, the two terms 

(S^Clw)* (H( j.w) ,.F(jw) ) ) 

N 2 (jw) 

iro^ 

may become much larger in magnitude than F(jjw ) . The inverse 
transform then no longer resembles f(x) 1 , but contains many 
noise-like variations. 

The difficulty seems to be that the form of the 
filter is fixed in the case .of inverse . filter. If* however, 
the filter is obtained by an estimation procedure where soma 
measure of the difference between the restored output and the 
original input is minimized in a statistical framework, some 
of the drawbacks of the inverse filter can be avoided;,. The 
most common measure is the mean square error between the input 
and the output. The corresponding filter is known as Wiener 
filter. 
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Here we will assume that the original input and 
both the additive and multiplicative noise components i.e* 
f(x), n 2 (x) and (x) respectively belong to zero -mean, 
mutually independent random processes. 

As described, by equation (2.1), the output c(x) 
is given by 

g(x> * / h(x-x f )f(x f )dx* +n 1 <x). S h(x-x‘ ) f (x‘ )dx* + n 2 (x) 

This equation can not be solved directly when noise 
is present, since both n^(x) and n^Cx) are not known to us 
though their statistical properties are assumed to be known. 
Instead, one can try to minimize the mean. eq. error 

a 2 « E { ( f(x) - f(x)> ! .. (2.5) 

in order to get the least-rq-estimate f(x) of f(x), given 
g(x) (E. denotes the expectation value). The unconstrained 
solution of this minimization, is evidently the conditional 
expectation of f(x) given g(x) (9 ), which is, in general, a 
non-linear function of g(x). Also one requires the joint 
probability density over the random processes f(x) and g(x) 

. I 

• A. 

making the process more complicated and inconvenient. 
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The problem becomes simplified if it is assumed 

A 

that the estimate f(x) is a linear function of the gray- 
levels in g(x). 

The corresponding estimate is known as the linear 
least square estimate. 

Mathematically, this is expressed by 

'fCx) = /m(x,x') g(x*> dx‘ . .. (2.6) 

where m(x,x‘) is the weight to be applied to g(x) at x =x‘ 

A 

for the conputation of f at X. We assume that both f (x), 
n-^Cx) and ^ (x) are second order stationary processes. In 
such a case, the weighting function depends only oa (x-x‘) 
and we have 

A 

f(x> = S m(x-x* ) g(x') dx‘ .. (2.7) 

end so 

* t 

e 2 * E t (f(x) - / m(x^x* )g(x* )dx‘ ) 2 J i. (2,8) 

2 

e can be minimized by several standard procedures (e.g. va- 
riational calculus techniques) and one obtains;, the following 
integral equation ( 4 ). ‘ 

Jm(x-x*) R gg (x , -7)dx’ = R fg <x-r) .. (2,8) 
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where R^(x-7) = S(g(x)g(7)> = Autocorrelation of g(x) . 

and Rfg(x-7)=» E(f(x).g(7)) = Cross-correlation of 

f(x) and g(x). 


From equation (2.8) , one finally obtains (2 ) the filter func- 
tion (in frequency domain) M(jw) as 


M(j,w) 


1 

|H{jw) 

I 2 35 f (w) 


™ H?jw) 

|h( jw)j 

2 $ f (w)+( |h(Jw) 

| 2 .$ f (w)>*a n (w)+a n (w) 
i 2 


.. (2.9) 


(*** denotes convolution.) 


where $£(w) * Power-spectral density of the input object f (x) , 

35^ (w) « Power-spectral density of the multiplicative 
noise n^Cx)*. 

and $ n (w) = Power-spectral density of the additive noise (x) . 

2 

(The PSD (Power spectral density), by definition, is the 
Fourier transform of the corresponding autocorrelation function). 


The M(jw) as given by equation (2.8) is the well-known Wiener 
filter for the model under discussion. 
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It is observed that the Wiener filter is indeed the 
filter multiplied by a noise smoothing factor M* (w) 

|H(jw)| 2 ffi f (w) 

Jh(jw)| 2 «B f (w) + ( |H(jw)j 2 <5 f (w))*« n (w)+C n (w) 

X 2 

.. ( 2 , 10 ) 

In noise-less environment i.e. when $ (w) = (5 (w) = 0, 

n l n 2 

the Wiener filter is exactly the inverse filter. A close 
examination of equation (2.9*) shows that the denominator 
actually represents the PSD of the recorded image and thus can 
be computed by any of the standard power-spectrum estimation 
techniques. However, the evaluation of the numerator requires 
the explicit too wl edge of a f (w) i.e* object PSD and thus one 
requires an object prototype for the computation of $£(w) . 
Mathematical convenience is the main reason for the choice of 
the minimization of m.s.e* (nreaa square error) as the crite- 
rion for filter design and unless a clear and reasonably 
accurate model of bumart visual process is established, this 
is likely to be a popular criterion. 

However, as already explained, the filter requires 

4 * 

an extensive knowledge of the PSD of the object, or its auto- 
correlation function* Though the PSD of the noise can be 


inverse 

where 

M« (w) 
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obtained from a suitable model, that of the object is essen- 
tially an intrinsic property of the object itself. In 
addition,. Wiener filter produces the optimum estimate only on 
an average. Since it solves an unconstrained minimization 
problem, false details due to negative pixel values are also 
encountered (2,7) . 

ITERATIVE WIENER FILTER (CONTINUOUS) : 

An iterative algorithm has been proposed and studied 
by Ramakrishna R.S. (2,7) which dispenses with the extensive 
knowledge of the object PSD for the computation of the Wiener 
filter (Fig. 2.2) . 

With an initial guess (w) of the object PSD, 
the Wiener filter is obtained as 


M(jw) « — — — 
H(jw) 


|H(jw)| 2 $ 1£ (w) 

-i m — p — m t m ■ mmm m*m mm • ym m mmmm mu ■ mrnmmm mm m 

|h ( j w) j 2 J5 lf (w) + ( jH(jw)| ffi lf (w))*<5 n (w)+35 n (w) 


|H(jw>| 2 .as lf (w) 
ffl lf (w) .H* ( jw) 


.. ( 2 . 11 ) 

1 

( ]H(jw)} 2 as lf (w))*® n (w) as n (w) 

• +— 2 : 

as lf (w).H*(jw) $ lf (w).H*(jw) 






\ 
i 
\ 

\ Stop j 
\ ' 


Fig. 2.2 Iterative algorithm 








If (w) and ? i (w) be the PSD-s corresponding to 
the restored and recorded images* then they are related by 

3Sq (w) « |M(jw)| 2 $ ± (w) . . (2,12) 

where ^(w) = }H(jw)f 2 <5 f (w) + ( |H(jw)| 2 S f (w>*£E n ^(w) 


+ 3L (w) 


n 


( 2 . 13 ) 


As is shown below* with an appropriate guess of the 
object PSD for the computation of the Wiener filter, the PSD 
of the restored output closely resembles that of the original 
object. 

From (2.11), (2,12) and (2,13), we get 


<E 0 (w) = 


|H(jw)j 2 a f (w)+( |H(jw)f 2 ai f (w))*(E n (w)+a n (w) 

_ 2 


t 

— 


jn(jw)j 2 as lf (w) ( jH(jw)j 2 $ lf (w))*$ n (w) a n (w) 


ffi lf (w)H*(jw) 


(w) ,H* (jw) 


1 + 


ll 


®lf (w) ,H* (jw) 


( |H(jw)j as f (w))*as n (w)+as n (w) 

1 + 1 2 

teH y iLgj f , ( a} 


l + 


( |H(jw)| 2 a lf (w))*as n (w)+ffi n (w) 

■ ' ~ 1 : , '' 3 

}H(jw)| 2 «i f (w) 


2B f (w) 
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G(w) *3i£(w ) 


.. ( 2 . 14 ) 


where , 


0(w) 


( jH(jw)j 2 s f (w))*ffi n (w>+ ® n (w) 

1 4. ^ % 


( ^(jw)| as lf (w))*a n (w) + « n (w) 

1 1 


n 2 


1 + 


.2 


|H(jw)j ffi lf (w) 


.. ( 2 . 15 ) 


If both |H(j,w)| 2 ^^(w) i.e. PSD of the degraded object and 

. 2 2 

]H(jw)J ffi lf (w) are reasonably large oonpared to ( ( jH(jw)j 

<t f (w))*(E n (w)+as n (w))and(( jH(jw)j (w))*5 n (w )+® n (w))res- 

12 12 
pectively at each w (at least over most of the frequencies 

under consideration) , then 0(w)— 1 and so a> o {w) Czi. 3i^(w) , 

(Note that C( jH(jw)l 2 $.g(w))*$ n (w) -HS^ (w)) actually represents 
the PSD of the component of the degraded output contaminated 
by noise) which shows that the PSD calculated from the first 
iterate of this procedure is very close to the object PSD 
provided the initial guess is reasonable. If the input SNR 
is known to be - quite high, one can coirpute the PSD of the 
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recorded object and use it as the first guess. Also one 
can start with an appropriate model of the input PSD and 
calculate the parameters associated from the recorded object 
itself (e.g. in many situations, it is quite reasonable to 
use a lawpass filter function as the object power— spectral 
density, or an exponential autocorrelation function) . In - 
further steps of iteration , the PSD is calculated from the 
iterates themselves. Proceeding in this manner, one can 
reach a stage when significant improvement is reached and in 
an interactive situation, one can step the process of itera- 
tion on obtaining reasonable restoration. This is well 
supported by corrputer results in the case of one-dimensional 
continuous signals in ( 2 ) and for discrete generalised 
Wiener filtering of two-dimensional data as described in 
the next chapter. 



CHAPTER 3 


• discrete ITERATIVE GENERALISED WISHER FILTER 

In this chapter* we will consider the restoration 
of discrete signals. The corresponding model of the imaging 
and recording system is explained below (Fig. 3*1) . 



Fig. 3.1 Discrete image formation and recording system. 

The following notations are used in this model. 

x '= A M-element data vector representing the original 
discrete data sequence. 

H = A M X M transformation matrix corresponding to the 
\ degradation process of the imaging system. For a 
space-invariant system* H is a circulant matrix. 

=s A M X M diagonal matrix* where the diagonal elements 
belong to the multiplicative noise process. 
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= A M-element additive noise vector, 
f = The M-element, degraded and noisy output vector. 

(Note that in our notation, lower-case letters are used to 
represent vectors and upper— case ones for matrices.) 

The output f is given by 

£ = Hx + Hx + n 2 •* (3.1) 

As in the case o£ continuous Wiener filter, here • 
also we assume that the input x and both the additive and 
multiplicative noise corrponents belong to zero -mean, mutually 
independent random processes. 

The discrete Wiener estimation of the input obj act 
essentially consists of finding a M X M filter matrix A which 
minimizes the mean square error between the original data 
sequence x and the restored sequence 'x. The mean square error 
e between x and x is given by 

e 2 - S ( (x-x) ^(x-x) ) .. (3.2) 

2 

A direct minimization of e as given by equation 
(3 *.2) leads to a filter matrix A which can be used to filter 
the recorded data in space-domain. However, it is possible 
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to perform the Wiener filter operation in different transform 
domains utilising unitary transforms and it is thereby possible 
to significantly reduce the coirputational complexities by 
selective confutation. This is known as ‘Generalised Wiener 
filtering* ( 1 ) . 



Pig, 3,2 Generalised Wiener filtering. 

Pig, 3,2 shows the corresponding block-diagram of 
the generalised Wiener filter. First# a M X M unitary trans- 
form T is applied to the recorded data vector f# followed by 
the multiplication by the M X M filter matrix (in the transform 
domain) A and finally# the restored output 'x is obtained by 
talcing the inverse transform of the filtered data. The restored 

A 

output x is given by 

£ - T”* 1 A T f .. (3,3) 

where T -1 » T* fc .. (3f.4) 

since T is a unitary transform. 


Let 


B = T -1 A T 


(3.5) 
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Nota that. B is related to A by a similarity- trans- 
formation (equation (3.5)) and so there is a one-to-one corres- 

2 

pondence between A and B. As a result, one can minimize e 
(equation (3.2)) w.r.t. B in order to compute the Wiener 
filter matrix A. Therefore 

A 

x « Bf .. (3.6) 

From equation (3.2) and (3.6), 
e 2 * E ( (x-x) ^(x-^) ) 

= ECxSc) - 2. E (x^x) + E (x%c) 

= E (x^*x) - 2E (x t Bf) .. (3.7) 

2 

Now B must be oho sen in such a way that e is minimuttr. We 

2 

differentiate e by the elements of the matrix B and equate 

2 

them, to zero i.e. we make v E (e ) * j where £0 j is the 
null matrix and ^ is the gradient operator, i.e. 
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It is to be noted that ^ is a linear operator 
and so can be applied before taking the expectation value. 

We have the following identities s 

V b (^B £> = xf* .. (3.8) 

57 B (f t B t Bf) =■ 2.B ff fc .. (3.9) 

and P^Cx fc x) *0 .. (3.10) 

Prora equations (3.7), (3.8), (3.9) and (3.10), 

V B (e 2 ) - 0 

or, e[V b (*"*)} - = 0 

which leads to 

E.(Bff fc ) = E(xf t ) 

OIT/ B*E (ff*) » E(xf t ) ... (3.11) 

Now from equation (3.1), 

f = Hx + Hx + n 2 

Hence E (xf^ = E. (xx fc ) ^ + E. (xxVl^) + E .(xa^) 

Now, E (xx**) «* R x l.e. autocorrelation matrix of the random 
process x. 
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E (x^2 *> -£0> since x and ^ are mutually independent 
zero-wean random processes. 

and E.CxxVi^) = £o7j, since x and are also mutually in- 
dependent, zero-mean random processes, 

(This is evident if one considers any element, say (i, j ) th 
element of E.Cxx^H^), since BCxxVl^)^ * E(E 


Thus E (xf^) « 

Also, E (ff^) = E ( (Hx+N^Hx-H^) (Hx-H^Hx-tt^) 


.. ( 3 . 12 ) 


as E (Hxx t ’H t ) + E (HxxV^) + E (Hx n^) 

+ E (N 1 Hxx t H t ) +E (^HxxV^) + ECN^Hx n^) 

+ E (n 2 x t H t ) + E. (n^H^) + EO^n^ 

As before E (ixx^) * E = 0 ; 

E(n 2 » 2 ^ * Rjj i.e, the autocorrelation matrix of H 2 » 

Also, since x, and n 2 are mutually independent, 
zero-mean random processes* the cross-terms E (Hxx^H^J^ ) , 
E(K 1 Hxx t H t ) , E(N X Hxn^ and E O^xV^) are zero (For details, 
see the Appendix ) . 
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Now, ECNjHx^H^) 


±S 


B * ( l H±i *k N 1 j:} H Jk 5 

* S *^ Ml ii Nl jj lk H il *1 *k H jk 5 


ii jj 

* Nl ij^ - ! a Hjx ^ H j* 


(since x and are mutually independent random processes,) 

® ij 


( (H R x H^) 


Rjj In the autocorrelation matrix of the multiplic- 


where 

ative noise process and ® denotes the Hadamard product of 
matrices , defined by (Ad B ) = a^j b^ , A and B being 

two xnatriees and a^j and b^ their (k,j)th elements respec- 
tively. 


Thus S (ff^) = H R X H^ + 0 R^, + R^ 

1 2 

From equation (3.5), (3.11), (3.12) and (3.13), 
B.E(ff t ) * E (xf fc ) 

or, B.(H + (HR^) 0 ^ + Rjj^) = R x H fc 

or, B - R x H t (H R x K fc + (HR Jt H t ) 0 + Rj, l" 1 

3. 2 


M (3.13) 
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which leads to 

A = T (HR x H t + Q ^ }‘’ 1 }T’* 1 .. (3.14) 

The matrix A, thus obtained from equation (3.14) is the 
generalised Wiener filter matr ix. Note that the restored 
output x is given by 

/V 1 

x = T A T f 

= T^T (R^ (HR x H t + (HR x H t ) O ~ 1 ) T^Tf 

*» (R x H t (HR^S (HR x H t ) 0 ^ +1^ )“V .. (3.15) 

i 2 

2 

which is independent of T and so the mean square error e = 

*t ’ 

E.(tx^oc) (x-x) ) does not depend on the transform chosen. 
However , the nature of the filter matrix A is strongly depen- 
dent oh T and thus one can choose a suitable transform in 
order to minimize the computational complexities in filter 
generation and filter operation! ( 1 ) • 

MINIMUM MEAN SQUARE ERROR : 

2 

The mean square error e is given by (equation 3.2) 
e 2 » E. ((x-xj^x-oc) ) « 
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or equivalently, 

e a Trace £ E | (x-x) (x-x)^j 

= Trace (Sta*)) - Trace (E .(*£*)) _ Trace (E&x*)) 

+ Trace (E (xx^ ) 

* Trace " 2. Trace (E(xx: t )) + Trace (E (xx^) ) 

(Since (*£*} is the transpose of C&t.) and so trace remains 
invariant. ) 

= Trace (R x ) - 2. Trace (E (Bfx t ) ) + Trace (ECBffV 1 )) 

(where S is given by equation (3.5)) 

* Trace tR x J ~ 2 «Trace (R^ fc (HR x #%HR x H t ) 0 + r^ )“ 1 hr^) 

1 * 2 * 

* Tra ce (R x H t (HR x H t + (HR^) e + R^ )~ 1 

1 2 

tHR x H + (HR x Ht ) ® *„ + \ > 

3 * 3 * 

(Hjyj* + (HR^H*) 0 Bj, + Hh ) _1 HR X ) 

3* 2 . 

(Using equations (3.12) and -(3.13). Note that 

(HR xH (HR^H^*) © r^ + r^ ) is a symmetric matrix) 

. . 3* 2 M 

= Trace (R x - + (HS x H t ) 0 ^ + R„ > _l hh*) 

1 2 

(3.16) 
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The minimum mean square error is given by equation 
(3*16) and as stated earlier, it is seen that the minimum m. s.e 
indep enaent of the transform chosen. Xn the absence of 
degradation and multiplicative noise i.e» when H = I (identity 


matrix) and R^ = £cfj , the minimum m. s.e. is given by 
e 2 = Trace (R* - !<*(“* + V* R x> 


= Trace j\(* x + V' j tR x + V 
= Trace (R^ + R,,) -1 R,,) 



.. ( 3 . 17 ) 


SUBOPTIMAL WIENER FILTERING ( 1 ) : 

The concepts of generalised Wiener filtering can 
be utilised to perform subop timal Wiener filtering in different 
transform domains. Here the problem is formulated as that 
of an unconstrained minimization of the m.s.e. with certain 
selected elements of the filter constrained to be zero. The 
main obj active is to achieve computational convenience with 
filter performance as close to the optimum filter level as 
possible. One selects the zero elements in such a manner 
that fast computation is possible. 

One example of this optimization is to use selected 
elements of the optimum Wiener filter matrix, given by equation 
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(3*14). For instance, only those terms which are of large 
magnitude can be retained. Another alternative will be to use 
only the diagonal and near diagonal terms of the optimum 
filter matrix. As can be seen from equation (3.14), the opti- 
mum. Wiener filter matrix A is given by 

A * T (R x H t (HR x H t +(HR x H t ) 0 + Rj. ) J )f 1 

= (T. (R x H t )T“ 1 ) (T( H R x H%( H R x H t ) © ) “ 1 T" 1 ) 

= (T(R x H t )T“ 1 )(T(HR x H t + (HR^ © + \ ) T" 1 )" 1 

1 2 

One can use a diagonal filter A', as a special case 
of this suboptimal filtering, where 

^ fT(R x H t )T ' 1 j(i,i) :| 

|^T(HR x H t + (HR x H t ) © R^ + JT** 1 )] U,i> 

12 i 

.. (3.13) 

In fact, when T is a real transform and the filter 
is constrained to be a diagonal one, the minimization of m.s.e., 
leads to the filter A* given by equation (3.13). This can 
be proved as follows. 

The restored output x, as before (equation (3.3)) 


is given by 



since T is a real unitary transform i„e. TT fc * i. x n equation 
(3.19), the filter A* is assumed to be a diagonal one i.e. 


A* * 


A‘(l,l> 0 0 

O A* (2,2) Q 


0 Q*. A* (M,M) 


f l 

j 


From equation (3*2) , 


e 2 * E -*■■£( x*oc)^(x~x)j 

= E (xSc) ■- 2 ECx^x) + E (x^x) 

* B (x\) - 2.E (xV A'T f) + E (fV A jt T T fc . 

A*T f) 

* E [x fc x] - aE^TXj^A'.CTf)! + E £(Tf) t . A^A*. (Tf)j 

.* ( 3 . 20 ) 

2 

The filter A* should be such that" e is minimized, 

\^,(e 2 ) * O, where in this particular case is given by 
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%« (S > - 


Let 


■^T^TTT 0 0 •• 


0 

0 


0 

0 


0 0 


X 1 - lb , 


'j (e 2 ) 
dA» <m,m) 


.. ( 3 . 21 ) 


*= Tf , 


.. ( 3 . 22 ) 


Then, from (3.20), (3.21) and (3.22), 

e 2 = E (x'ix) - 2E (x^A*^) A ft A*^) 

Now * XL { s <*>*4 ■ ft , | EC J t * i *' . 

(where, for a vector R, denotes the ith element ) 

“ D ( S <*!*!*>] 

£ where D is an operator, which,, when operated on a matrix 
takes the diagonal elements and produces the corresponding 
diagonal matrix. J 
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•Dfl.R^Vj ... (3 . 

Similarly, 


23 ) ( Prom equation 3.12) 

• - ' - t , , t„ 


^ i s < v a‘X>] 

I s ( J x A,(1 ' i)2 *t M 


i J 


2.D £ E (A* f 1 £ L t ) J 
2-.D ( E(A* T fft T^J 


** 2; • I> [ A, T. E (££*) T^J 

= 2.D [A' I (HR x H t + (HRH*) 0 ^ + „ , ^ 

(From, equation ( 3 . 13 )) 

= 2AM) {KHIy^ + (HR^) 0 1^ + ^ ) T tj 

.. (3.24) 

(Since A* itself is a diagonal matrix) 

Hence from equations (3.20), (3.23) and (3.24) 

we get 

V A , < e2 ) * o 
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or. 2.A-.D { T (HR x H fc + (HR^) © ) t* j 

1 2 

= 2,d[t. fyiVj 

which leads to the following equation for the diac 
marts of the diagonal matrix A* s 

A* (1,1) = 

£ Since T“^ = | 

which is the suboptimal filter already described by equation 
(3.18) . The resulting filter is computationally convenient 
because (i) it is a diagonal filter and so filter matrix 
multiplication is less complicated, (ii) it avoids matrix 
inversion. 

Iterative generalised Wiener filtering : - The discrete, 
generalised Wiener filter, like the continuous one can be 
criticised on the following grounds * 

(i) The filter, as given by equation (3.14) requires 

an extensive knowledge of the autocorrelation 
matrix of the original object. . Though the autoco- 
rrelation properties of the noise can be obtained 




from a suitable model, that of the object is 
essentially an intrinsic property of the object . 
itself. 

(ii) The Wiener filter produces the optimum estimate 
only on an average and the® is no guarantee that 
every filter operation will produce the optimum, 
restoration. 

(iii) Since the Wiener filter solves an unconstrained 
minimization problem, false details due to negative 
pixel values are also encountered. 

The same iterative technique, which was proposed 
to avoid the above mentioned problems in the case of conti- 
nuous Wiener filter, works out here also. 



Fig, 3.3 ,/ Iterative generalised Wiener filtering 
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Thus one starts with an initial guess R of R v * 

x o x 

obtains an estimate* uses it to compute a better (hopefully) 
obj ect autocorrelation matrix and proceeds in an iterative 
manner till significant improvement results. Also note that 
the iterative procedure can be applied in the case of 
suboptimal filtering also* since* from equation ( 3 . 18 )* it is 
seen that the diagonal filter A* also requires the precise 
knowledge of R^ 

The iterative Wiener filter does not require the 
extensive information about object autocorrelation matrix. 

In addition* non-negativity of the pixel values may be 
introduced in the steps of iteration. Moreover* the whole 
algorithm is tied to the given recorded image and so 
reasonable restoration. may be expected for each object (and 
not just on the average). 



CHAPTER 4 


SIMULATION. & RESULTS 

SIMULATION OP THE DEGRADED AND NOISY IMAGE s 

The proposed scheme for iterative, generalised 
discrete Wiener filtering was implemented- with two 16 x 16 
images, blurred by two row and column convolution matrices 
H_ and H where 


H C t=H R 


If P is the original picture matrix, then the 
blurred matrix P‘ is given by. 

p« = Hj^.P «H C •• 

If v be the vector obtained from the picture matrix 
P by stacking the elements columnwise i»e«, 
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V - |p (1/1) / P(l,2), .... p (1>M) t P(2,l) , P (2,2) , 

.... P(2,M), .... P(M,M-1), P(M,M)J .. (4.2) 


than the corresponding degraded output vector v* (which 
also be obtained from; P * by similar process of columnwise 
stacking) is given by (3 ) 


v* =* H.v 


.. (4.3) 


where. 


n = h r a h c 


(4.4) 


St denotes left direct (Kronecker) product of two matrices 
i.e. if A and B be two MXM matrices, then 


a a b * 


B u[ A ] B 1 2 M 


B 


lM 


W 


B 2 l[ A ] 


B 22l A ] 


• • " B 2mW 


b mi[ a ] 


B M2[ A 1 


B. [A] 


MMi 


If A and B are circulant matrices/ then A SB is a block 
circulant matrix. Note that, in general in an imaging system. 
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one may not be able to separate out the row and columnwise 
degradation processes. In such a case H (equation (4*4)) 
can not be expressed as a direct product of row and columnwise 
^s^t’ndation matrices and the blurring is described by equation. 
(4.3). Thus equation (4.3). is the most general representation 
of image degradation. 

Zero-mean, unit-variance, white, Gaussian noise is 
added to the degraded picture. Moreover, the picture is also 
corrupted by zero-mean, white, Gaussian multiplicative noise. 
The noise auto-correlation matrix is clearly a diagonal matrix, 
since for a zero-mean, white noise process, 

E(n ± nj) = •• C 4 * 5 ) 

where n^ and n^ are the i-th and j-th sanples of the white 
noise and 0 is the noise variance. 

OBJECT AUTOCORRELATION! AND ITS COMPUTATION s 

The autocorrelation R(i,k) of a discrete random 
process x(i) is given by 

R(i,k) * E(x(i).x(i+k)) .. (4.6) 

If the discrete sequence vector x, of which x(i) is the i-th 
conponent belongs to a second order stationary process 
£we define x £ (x(l), x(2), .... x(M) ) J , then 
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autocorrelation becomes a function of k. For such a stationary 
process, equation (4.6) modifies to, 

R(k) ■ E(x(i).x(i+k) ) .. (4.7) 

for any i. 

The object auto-correlation matrix R x is given by 


R x * E(xx t ) 



S (x(i) .x(l) ) 

E(x(i) .x(2) ) .... 

E (x(l)x(M) ) 

E(x(2) ,x(l) ) 

E(x(2).x(2)) .... 

• 

B(x(2).x(M)> 

• 

• 

E<x(H) »x(l ) ) 

• 

♦ 

E(x(M) ,x(2) ), . . . 

* 

1 

E(x(M)*x(M) ) 


Since E(x(i) ,x(j) )= E (x(j ) .x(i», R x is a symmetric 
matrix. From this and the assumption of stationarity of x, 
it is evident that R v takes the - following structure. 
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(ay the arrows, we mean that the elements R(0), R(l), R(2) 

^ (M— 1 ) are continued in the corresponding lines.) . 

Such a matrix is a symmetric Toeplitz matrix. For 
a two-dimensional object, both nowise and columnwise auto- 
correlations take part in the object . ottocor relation . With 
v being the vector obtained from the picture matrix P by 
stacking the elements columnwise (equation (4.2)), the 
corresponding autocorrelation matrix R w = E (w fc ) takes the 
following form. 



where denotes the i-th MXM block. Such.* matrix is known 

as a block-Toeplitz matrix. Here each block is a MXM Toeplitz 
matrix because of rowwise' stationarity and the blocks are 
arranged in Toeplitz form because of columnwise stationarity* 
Note that each block , i — 0,1, . •» (M-l), incorporates 

the autocorrelation between two rows at a gap of i. 
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Also note that if the object autocorrelation can 
be separated both row and colwnrwia* # then can be 
expressed as 

R w = R S a R 0 .. (4.10) 

where, as before & denotes the left direct product of matrices 
and R c and 
matrices respectively. 

The computation of the obj act autocorrelation matrix 
is based on the assumption of ergodicity of the data sequence. 
The sample autocorrelation function R(k) is computed by the 
following equation : 


Rp are the columnwise and rowwise autocorrelation 


R(k) 


1 

n 


n-1 

Z 

i=0 


x i x i+k 


.. (4.11) 


where n = M-k •• (4.12) 

(M is the length of the sequence.) 

Now# consider a two-dimensional sequence x, given by 
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x o,o 

*0,1 

• * • • X 0,M-1 

*1,0 

*1-1 

* * * # *1,M-1 


*M-1,0 *M-1,1 ‘ 


The corresponding autocorrelation R^,}^) i.e. 
autocorrelation between saitples at a separation (k^jk^) is 
computed from the following equation s 


RO^ # k 2 ) = 


1 

mn 


n-1 m-1 

if o j fo Xl/ J * Xl+k l * ' j +3c 2 


.. (4.13) 


where 



.. (4.14) 


n * M-k^ .. (4.15) 

It is to be noted that for stationary data, the 
computation of the autocorrelation matrix requires the 
evaluation of any row# or, column only, because of the sy- 
mmetric Toeplitz (or block-Toeplitz in case of two-dimensional 
signals) structure of the autocorrelation matrix. Thus one 
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can evaluate the autocorrelation values for different values 
of k (or k^* for. 2D data) using equation (4.11) (or (4.13) 
for 2D data.) and thereby obtain the autocorrelation matrix. 

The proposed scheme of iterative* generalised* 
Wiener filtering was applied to two 16 x 16 images. The 
initial guess for object autocorrelation was derived from the 
recorded data itself* With input SNR reasonably high and 
degradation reasonably less* the choice is quite appropriate. 
Also one can start with a suitable model for the object 
autocorrelation and compute the several associated parameters, 
from the given data (e.g. in many situations* it is quite 
appropriate to use a low pass filter function 

as the object PSD* or equivalently* an exponential autocorre- 
lation function). The mean of the object was computed from 
the iterates in each step of iteration and substracted from 
the iterates themselves* thereby generating a zero-mean input. 

RESULTS 

The simulation results corresponding to the two 
images are described below* 

(A) The corresponding image is shown in Fig* 4.1(a). 

The noisy and degraded image is shown in Fig. 4.1(b) . 
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Variance of additive, white, Gaussian noise = 1.0 
Variance of multiplicative, white, Gaussian noise = 0.1 
Initial mean square error = 18.47 
(i) Optimal space domain filter 


A * 



NO. Of 
Iterations 




Mean square 22.15 15.61 12.00 10.13 8.01 7.91 

error 


(ii) Suboptimal filtering 

P^i -1 ) u,D 

" 0 \ + I_1 ] 
1-1 t 

where T = T* and A is a diagonal matrix, 
(a) T = I 

No. of iterations 
Mean square error 


12 3 4 

11.01 10.00 9.89 9.76 
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(b) T *e Hadamard transform 


No- of iterations 

Mean square error 

1 

11.79 

2 

10.62 

3 

9.82 

4 

8.93 

5 

8.89 

(c) T « FFT 






No. of iterations 

1 

2 

3 

4 

5 6 

Mean square error 

13.26 

11.85 

11.03 

10.69 

10 .£>1 9.46 


The pictures corresponding to the different steps 
of iteration are shown in Figs. 4.2 - 4.5. 
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B) The original image and the corresponding noisy 

and degraded picture are shown in Pig. 4 . 6 (a) and Pig. 
4.6(b) respectively.. 

Variance of additive, white, Gaussian noise = 1.00 

Variance of multiplicative, white, Gaussian noise = o.lS 
Initial mean square error » 31 34 

* r x 


(!) Optimal space domain filter 



No. of iterations 1 2 3 4 5 6 

Mean square error 25.16 12.50 11 . 99 10.26 9.87 8.71 


(ii) Suboptimal filtering 

(T R T -1 1 (i,i) . . 

A(i,i) * ~ £— ~ — — -‘T’-xr—- 

[T(HP^ + (Hiy-T) €> Rj, + .*$ ) * J (i,±> 

—1 t 

where T * T* and A is a diagonal matrix. 

(a) T m I 


1 2 3 4 5 

16.31 15.71 15„53 15.17 15.15 


No. of iterations 
Mean square error 
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(b) T » Hadamard transform 


No. of iterations 

Mean square error 

1 2 3 

14.86 12.98 10.75 

4 

9.78 

5 

9.81 

(c) T » FPT 

No., of iterations 

1 2 3 4 

5 

6 

Mean square error 

16,45 14.86 13.52 13.57 

12.33 

10.37 


The corresponding pictures are shown in Pigs. 4.7- 


4.10 
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Hg» 4 o 6 (a) Original image s 

(b) No lay sad degraded image* 
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CHAPTER 5 


CONCLUSION 

An attempt has been made in this thesis to remedy 
some of the shortcomings of conventional Wiener filter through 
an iterative algorithm. The limited simulation studies of 
the algorithm for the case of generalised Wiener filtering 
have shown the effectiveness of the iteration strategy. 
Mathematical proof of the convergence of the algorithm is 
a formidable task and it is felt that more work is required 
to solve the problem. 

The restoration problem can be simplified with 
the aid of circular ts for space-invariant degrading systems. 
The concept of asymptotic equivalence of circulsnts and 
Toeplitz matrices can be used to approximate Toeplitz (Block- 
Toeplitz) correlation matrices by circulants (Block-circulants 
The latter permit the use of FFT-algorithms and the problem 
becomes computationally convenient. 

It is observed that for suboptimal filtering, the 
iterative algorithm works better in the case of processing in 
Hadamard domain than in FFT domain. This can be attributed 
to the fact that the pictures used for simulation are binary 
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in nature i.e. they have only two levels, high and low and 
so their representation and processing in Hadamard 'domain 
leads to less loss of information in comparison to PFT. 

It is expected that for pictures with continuous variation 
of gray levels, FFT-domain processing would work better. 
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APPENDIX 


In chapter 3 , we mentioned that 

(i) ' E (H xx? Nj) * [o] 

(ii> ECNjH xx? H 1 ') * £oJ 

. (iii) EC^ Hx n^) *£'03 

(iv) E(n 2 x t N^) = £oJ 

where, as described in chapter 3 , the original object x and 
the additive and multiplicative noise corrponents i.e. n 2 and 
N-^ are mutually independent, zero -mean random^ processes. Here 
we give a formal proof of these results. 

<i) E (Hxx^ N^. . 


" E ( l l *il X 1 *k N l j:j V 

* ( n H ii- E tx x x^), ECNj^ ) Hj ) 


= 0 , since the object x and both and N 2 


zero-mean, mutually independent random processes. 
Hence ECHxxV?^) * [OJ 



(ii) 


(iii) 


Civ) 


E (N, H X3C H t ) . . 


=E{ ! t Ni ii x i ** V 
= ik EtNl ii )- H a- E(x i.v v 


* 0 r Hence E (N^H xx^ H^) = £o ] 
E(N X Hx 

■*' k^ii H iK- ^”2,’ 

* C ? E(N i 11 ) - H ix* E <*k>- «<»a >> 


Hence E (N 1 Hxn 2 t ) = [qJ 


E (n 2 x t H t N 1 ) 


=* E(S n 


, 2 *|f' N i H., ) 

k 1 ^ 1 j,j Jjk 


(S ECn^)^ E^). Eft^ )., H Jk ) 


* o 


Hence ECn^H^) » £oj 





